Synchronisation of stochastic oscillators in biochemical systems 
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A formalism is developed which describes the extent to which stochastic oscillations in biochemical 
models are synchronised. It is based on the calculation of the complex coherence function within 
the linear noise approximation. The method is illustrated on a simple example and then applied 
to study the synchronisation of chemical concentrations in social amoeba. The degree to which 
variation of rate constants in different cells and the volume of the cells affects synchronisation of the 
oscillations is explored, and the phase lag calculated. In all cases the analytical results are shown 
to be in good agreement with those obtained through numerical simulations. 
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I. INTRODUCTION 

Oscillatory behaviour is observed in a great variety 
of biochemical systems, over a wide range of time pe- 
riods [TJ These oscillations have been modelled ex- 
tensively, using both deterministic and stochastic frame- 
works |3H6]. Most of this modelling is carried out at the 
single cell level. However, many processes display coher- 
ent behaviour over a population of cells. This requires 
that the oscillations in the individual cells must influence 
each other. One reason for this is that random fluctua- 
tions, intrinsic to these systems, can introduce random 
phase lags. This means that a population of isolated 
cells demonstrating oscillatory dynamics will not remain 
synchronised over time, even if they are synchronised ini- 
tially. Therefore, for coherent, collective behaviour (such 
as intercellular signalling), some form of communication 
between the cells is necessary. 

Before proceeding further we need to discuss what is 
meant by the term 'synchronisation'. This is necessary to 
do this since its precise definition varies. For the stochas- 
tic time series of the kind that will interest us here, the 
approach of Pikovsky, Rosemblum & Kurths will be the 
most relevant. They define synchronisation as "an ad- 
justment of rhythms of oscillating objects due to their 
weak interaction" [7] . Using reaction systems as exam- 
ples, we will show the forms that the interaction can take, 
and what 'weak' means in this context. 

One place in which the role and mechanisms of syn- 
chronisation have been systematically analysed is the 
study of the dynamics of glycolytic metabolism. Richard 
et. al. have performed experiments in yeast cells, and 
observed oscillations both in individual cells and pop- 
ulations of cells [H [9]. The authors conclude that the 
coupling of the cells via an extracellular metabolite (pro- 
posed to be acetaldehyde) causes the oscillations of the 
individual cells to synchronise. In j^, two cell popula- 
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tions, originally oscillating 180° out of phase, became 
synchronised upon coupling, and a time to synchronise 
was measured. 

Several mathematical models have attempted to cap- 
ture and explain this behaviour. One of these is due 
to Wolf & Heinrich, devised in [10] and refined in II Ij. 
They begin by describing a deterministic, one-cell model, 
where the cell is connected to an extra-cellular compart- 
ment, which exhibits either a stable fixed point or a stable 
limit cycle, depending on the values of the reaction pa- 
rameters. Wolf & Heinrich then extended the model to n- 
cells, where the cells are coupled by exchanging molecules 
with a shared extra-cellular compartment. In [TT], the 
authors show that in a two-cell model, synchronisation 
with either zero-phase or non-zero phase lag is observed, 
depending on the reaction parameters chosen. Follow- 
ing the experiments by Richard et. al., the authors also 
mixed two populations of cells which, before mixing, were 
internally synchronised but out of phase with the other 
population. Wolf & Heinrich found that, once mixed, the 
two populations did indeed synchronise, although this 
took considerably longer to achieve than it did in the 
experiments. This may suggest that the form of the cou- 
pling, or the chemical species chosen to be responsible 
for the coupling, may not be the correct one. In many 
models it has been found that the form of the coupling 
chosen in the model can significantly alter the dynamics 
observed, as found in [T^ I13j . 

Another system in which synchronisation is often stud- 
ied involves oscillations in the concentration of adenosine 
3', 5'-cyclic monophosphate {cAMP) in the amoeba Dic- 
tyostelium discoideum (2 1^ . These social amoeba feed 
on bacteria in the forest soil. During the onset of star- 
vation, the cells alter their behaviour to aggregate, and 
become able to relay cAMP signals, in the form of oscilla- 
tions. In this way, populations of the Dictyostelium cells 
form large aggregates. This process culminates with the 
formation of a fruiting body which can disperse spores, 
leading back to the start of the life cycle . One reason 
why this is much studied is that many of the compo- 
nents involved in the cAMP signalling have correspond- 
ing components in mammalian pathways, hence the de- 
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sire to understand this system in more detail [T3] . 

In this paper we will present an approach to the in- 
vestigation of synchronisation in models of biochemical 
processes which is analytic and where the models are 
stochastic in nature. The vast majority of studies of such 
systems found in the literature are described by systems 
of ordinary differential equations (ODEs), which show 
limit cycle behaviour. However it is now well established 
that stochastic effects will be significant even in sys- 
tems containing a large number of molecules [SJ Uni dZ] . 
Stochastic oscillations will not only occur at a single 
frequency, as for limit cycles, and this will have conse- 
quences when investigating synchronisation of these os- 
cillations. What studies there have been of synchronisa- 
tion of stochastic oscillations have been purely numeri- 
cal [18j , and have not considered the possibility of phase 
lags between oscillations in different cells. However, the 
key aspects of synchronisation in stochastic biochemical 
models become particularly clear when an systematic an- 
alytical treatment is given. 

The outline of the paper is as follows. In Section [iT] we 
will introduce the formalism we shall be using and give a 
simple example to illustrate its use. A model of cAMP 
oscillations in Dictyostelium discoideum will be analysed 



in Section HI to show the application of our methods 
to a more realistic, multi-cell model. In particular, we 
show under which conditions a phase lag between oscil- 
lations in neighbouring cells is introduced. We end with 
a summary and suggestions for further work. 



II. METHODOLOGY 

In this section we will describe the formulation of the 
process in terms of a stochastic model, and then go on to 
describe the methods we use to see if synchronisation oc- 
curs. Fortunately, the quantities we need to calculate to 
investigate the question of synchronisation are the same 
as those required to study stochastic oscillations, and so 
in many cases much of the calculation will already have 
been carried out. 



A. Formalism 

The reaction system is described by a chemical master 
equation, which is written in terms of the integer pop- 
ulations of the chemical species present in the model, 
which we write as n = (rii, n2, tik), if there are K 
species present in the model. This describes the time 
evolution of the probability for the system to be in state 
n at time t, which is written as Pn{t). Transitions from 
one state to another are caused by the chemical reac- 
tions. We write T{n'\n) as the transition rate from state 
n to state n' = n — where i/^ = (i^i^, . . . , vj^^) is the 
stoichiometric vector corresponding to reaction /i. The 



master equation then has the general form [121 120] 



(1) 

-T^{n + u^\n)Pr,{t) 



To make progress analytically, an approximation 
scheme must be employed. The crudest scheme is to ig- 
nore fluctuations completely, and simply write down an 
equation for the average 



,{t)) =^n,Pr,{t), i^l,...,K. 



(2) 



In the limit where both the number of molecules in the 
system and its volume, V, becomes large, we describe 
the system in terms of the concentrations of the chemical 
species, Xi — limv^ooini) /V . From Eqs. ([I]) and ([2| one 
finds the macroscopic, deterministic, equation for Xi{t): 



A,ix), Aiix) 



dt 

where ff^{x) is defined by [2Ql 



(3) 



f,{x)^ Jim T,{{n}+,.,\{n})/V. (4) 

It is straightforward to calculate the function Ai{x) from 
the chemical reactions and the reaction rates, and so we 
can determine whether at large times the system tends to 
a fixed point, a limit cycle, or some other, more complex, 
type of behaviour. 

However, as mentioned already, oscillations seen in bio- 
chemical reactions may be stochastic in nature, and so 
not be described by the deterministic equation ([s]). To 
study these we cannot completely ignore the dynamics 
of the fluctuations contained in Eq. ([I]). The simplest 
way of doing this is to simply include stochastic effects 
(or 'noise') as linear deviations from the deterministic re- 
sult. This is the linear noise approximation (LNA) and is 
especially simple to implement if the system has reached 
a stationary state which corresponds to a fixed point 
of the deterministic equation (|3]). The LNA then con- 
sists of working with the variables Xi{t), writing them as 
Xi — X* + V^^l'^^i. The LNA assumes that Y is large, so 
we keep only linear terms in when they are substituted 
into the master equation ([T]). Here the asterisk signifies a 
fixed point of Eq. ([s]) and the V~^l'^ reflects the Gaussian 
nature of the approximation |19j . The equation satisfied 
by the linearised fluctuations ^ is found to be (for details, 
we refer the reader to the literature |51 [TBI EOj ) 



d^ 
dt 



K 

^M,,0+?7.W, i = l,...,K, (5) 



where rji it) is a Gaussian noise term with zero mean and 
with correlator {rji{t)rij{t')) = Bijd{t — t'). The matrix 
M is the Jacobian at the fixed point and the matrix B 
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can be calculated in the same way that the function A 
was: 



M,. 



Bi 



E 
E 



(6) 



Since Eq. ([5| describes fluctuations about the station- 
ary state and is linear, we can analyse the equation us- 
ing Fourier transforms. Denoting the temporal Fourier 
transform of as we may write Eq. ([5| as 

^^{^) = ^ij^{^)V3i^) where ^ij = -Mij-iuiSij. This 
then allows us to calculate the power spectral density ma- 
trix (PSDM) 



P,,M^(|.(a;)|*M) 



K K 

EE 

;=i m=l 



(7) 

In earlier work on analysing the nature of the stochastic 
oscillations about a fixed point, only the diagonal entries 
of this matrix (the 'power spectra' Pa) were of interest 
OEl]. Here we shall also be interested in the off-diagonal 
entries, but as mentioned in the Introduction, if Eq. ([t]) 
has been calculated so as to obtain the power spectrum, 
the off-diagonal elements can be immediately read-off. 

An application quite similar to the one we are dis- 
cussing here is that described by Rozhnova et. al. who 
used the off-diagonal elements to show how stochastic 
oscillations can become synchronised [22]. They con- 
structed an epidemiological model of a network of cities, 
where a disease can be transported between cities due to 
movement of infected commuters. We will use a similar 
approach but instead of cities and people, we will work 
with cells and molecules. Unlike the diagonal elements 
— the power spectra — the off-diagonal elements of the 
PSDM will in general be complex. Often, these elements 
are normalised, by using the relevant power spectra. This 
quantity is then known as the complex coherence func- 
tion (CCF) 



(8) 



As this is a complex function, it can be expressed in 
terms of a magnitude and a phase. 



\C^M\ 



(9) 



<j)ij{uj) = arctan 



Re(C,,(u;))^ 



arctan 



lm{P^^ 
Re(P,,M)^ 



(10) 

The magnitude of the CCF tells us the similarity between 
two signals, as a function of oj. The phase of the CCF 
tells us the phase lag between two signals as a function of 
the frequency oj [24]. For example, if the two signals are 
in phase the CCF will be a real function. To clarify the 
use of these concepts we apply them to a simple example. 



B. Simple example 

The biochemical reaction scheme we will use to il- 
lustrate the methodology is a toy system, but many 
more complex reaction schemes contain molecular species 
which behave in this way. The system consists of two 
species A and _B, which may both lose molecules through 
diffusion out of the cell and also gain them through diffu- 
sion into the cell. Species A requires molecules of species 
B to sustain its numbers, and they become depleted when 
too many B molecules have been used up. This is essen- 
tially a predator-prey dynamics, and we will use this lan- 
guage to discuss it, since it makes the results more intu- 
itive. The sustained stochastic oscillations have been pre- 
viously investigated using the formalism described in Sec- 
21]. The actions of these species is described by 
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tion 

the following interactions and their corresponding tran- 
sition rates: 



B + E ^ B + B, 



T{ni,n2 + l\ni,n2) = bn2 



{N - ni - n2) 
N '■ 



A — > E, T{ni - I,n2\ni,n2) = dni, 
A + B ^ A + A, T(ni + l,n2-l|ni,n2) -Pi^^, 
A + B^A + E, r(ni,n2- l|ni,n2) =P2™. 

(11) 

The variable E is introduced to give the system a maxi- 
mum carrying capacity, which is the 'size of the system' 
N. That is, the total number of individuals at any time 
is fixed, ni + n2 + ue = N, where ni,n2 and ue are, 
respectively, the numbers of A, B and E molecules in the 
cell. Because of this we will use N as the large parame- 
ter employed in the LNA in this example. We have also 
only allowed for species A to diffuse out of the cell, and 
for species B to diffuse in (subject to existing concentra- 
tion). The relation He = N — rii — n2 has been used to 
simplify the transition rates above. 

Applying the methodology described in Section |II A[ 
the deterministic equations ^ are 



dxi 

dX2 



■ P1X1X2 — dxi, 

bX2{l ~ Xi - X2) - {pi + P2)X1X2. (12) 



To go beyond this, we use the LNA which is fully char- 
acterised by the two matrices AI and B, which can be 
calculated from Eq. Q , and which are given explicitly in 

The macroscopic equations have a single nontrivial sta- 
ble fixed point (or 'steady state'), and for the reaction 
parameters chosen here the eigenvalues of the Jacobian, 
evaluated at the fixed point, have an imaginary compo- 
nent and sustained stochastic oscillations can be observed 
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0.05 



0.10 



0.15 



0.20 



(w) = tan 



/3 + 7a;2 



(13) 



for this two species model. Here a, j3 and 7 are calculable 
combinations of the parameters 6, c?, pi and p2 defining 
the system. The phase lag is shown in Figure|3] along with 
simulation results. Figure[l]indicates that the oscillations 
are most significant over the frequency range 0.05 < uj < 
0.11. The phase tends to tt as w — > and w — > 00. 
However within the range 0.05 < uj < 0.10 the phase lag 
remains fairly constant, at around 2 radians. 

We end this section by remarking that we do not use 
the word synchronisation to describe the phenomena ob- 
served in this section. We will only use this term for the 
interaction between self-sustained oscillators. We do not 
have two self-sustained oscillators here, as neither species 
would continue to oscillate if the other were removed or 
held fixed. In the next section, we will look at a more 
realistic, multi-cellular model, where self-sustained oscil- 
lators in neighbouring cells become synchronised due to 
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FIG. 1. Power spectra of the oscillations observed in the two- 
species example calculated using the LNA. The red (solid) 
curve is for the prey and the blue (dashed) curve for the 
predator. The parameter values used are: b = 0.1, d = 0.1, 
pi = 0.25, p2 = 0.05, and N = 3200. 



around the fixed point |21J. This is illustrated in Figure[T] 
which shows the power spectra of the oscillations, calcu- 
lated using Eq. ([t]) for a particular choice of the reaction 
parameters. The large peak shows the existence of ampli- 
fied stochastic oscillations for a reasonably narrow range 
of values about a characteristic frequency [5T]. 

We use the complex coherence function (CCF), intro- 
duced in the previous section, to study the relation be- 
tween these two oscillators. In general, the CCF will be a 
complex- valued function. Figure [2] shows the magnitude 
of the CCF, showing a strong shared signal in the two 
oscillators. The analytical result obtained from Eq. (|9]), 
is compared with results from numerical simulations, ob- 
tained using the Gillespie algorithm [2F. The phase lag, 
defined in Eq. (10 1, is found to have the simple form 
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FIG. 2. Magnitude of the CCF Ci2{oj) in the predator prey 
model. The orange curve is the theoretical prediction and the 
blue dots are the results obtained from numerical simulation. 
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FIG. 3. The phase spectra in the predator prey model, show- 
ing the phase lag (in radians) between the oscillators as a 
function of lj. The red line is the theoretical prediction and 
the blue dots are the results obtained from numerical simula- 
tion. 



the coupling between the cells. This is what we mean by 
the 'weak interaction' between oscillators. 



III. COLLECTIVE BEHAVIOUR IN 
DICTYOSTELIUM DISCOIDEUM 

The reaction system we will now discuss is a model of 
oscillations in the concentration of adenosine 3', 5 '-cyclic 
monophosphate {cAMP), which are observed within the 
amoeba Dictyostelium discoideum. The model is due to 
Kim et. al. and in their work was used to study the 
synchrony of stochastic oscillations in the concentration 
of cAMP across multiple cells via numerical simulation. 



OJ 



5 



Previously, the pathway that produces these oscillations 
had been modelled using a set of non-linear ODEs, which 
exhibit limit cycle behaviour [5]. Starting from such a 
model, described in [23" , Kim et. al. [15^ show that con- 
sidering stochastic oscillations, in addition to those in 
the limit cycle regime, increases the robustness of the 
oscillations. That is, oscillations, whether deterministic 
or stochastic in origin, are observed over a larger vol- 
ume of parameter space than would be the case if only 
deterministic oscillations were considered. Kim et. al. 
also demonstrate that, in the deterministic model, the 
oscillations cease to be observed if the reaction param- 
eter values are slightly altered. The robustness of the 
oscillations in the model is desirable, as any parameter 
values present in a mathematical model of a biochemical 
reaction system will have an associated uncertainty. In 
addition to this, because the physical system involves an 
aggregation of cells, conditions {e.g. temperature) may 
vary from cell to cell, so the reaction parameters may not 
be identical inside each cell. Because of this, Kim et. al. 
made small perturbations to the reaction parameters in 
different cells to see if oscillations were still observed, and 
whether they still synchronised. 

In this section we will use the theoretical framework, 
outlined in Section 2, to obtain analytical results for the 
system in terms of the CCF. This provides information 
on the strength of coupling between the oscillations in 
neighbouring cells (through the magnitude of the CCF), 



as well as the phase lag between the oscillations, if any. 
The possibility of a phase lag, introduced due to the het- 
erogeneity of the reaction parameters in different cells 
was not considered in [TSj. These calculations will be 
performed whilst parameters in the model, such as re- 
action rates and cell volumes, are varied, to identify the 
impact of changing these parameters. To study these 
multi-cell models in a stochastic framework we will use 
some ideas from previous work, where the LNA was ap- 
plied to multi-compartment biochemical models |27j . 

Although we will be interested in a multi-cell model, 
we will begin with a one-cell model, which consists 
of a cell and an extra-cellular environment. The 
cAMP within the cell is denoted by cAAIPi, the 
cAMP in the extra-cellular environment is denoted 
by cAMPe. The other components present in the 
model are adenylyl cylase (AC A), protein kinase (PKA), 
mitogen-activated protein kinase (ERK2), the cAMP 
phosphodiesterase (RegA) and the ligand-bound cell 
receptor (CARl). The coupling between the two 
compartments is caused by the external cAMP binding 
to the cell receptor. To label the species we use: 
(ACA, PKA, ERK2, RegA, cAMPi, cAMPe, CARl) 
(xi, X2, xa, a;4, xs, xg, xy). The reaction scheme, together 
with the reaction kinetics and numerical values of the 
reaction rates are given in Table 1, which can be used to 
calculate the matrices M and B using Eq. 



Reaction Scheme 


Reaction 






Rate Constant Value 


CARl ACA + CARl 


kiXr 


1^11 = +1 


ki = 1.96min~^ 


ACA + PKA PKA 


k2XiX2 


1^12= -1 


k2 = 0.882/iM-imin-i 


cAMPi PKA + cAMPi 


k-iX5 


V2-i = +1 


fca = 2.55min~^ 


PKA-^^ 


k4X2 


V24 = —1 


k4 = 1.53min~^ 


CARl ERK2 + CARl 


ksxr 


1^3 5 = +1 


fcs = 0.588min"^ 


PKA + ERK2 PKA 


kQX2Xz 


V-i(, = -1 


fee = 0.816/iM"^niin"^ 


RegA 


kj 


U47 = +1 


k7 = 1.02^iMmin"^ 


ERK2 + RegA ERK2 


ksXsX4 


U4S = —1 


kg = 1.274AiM-^min"^ 


ACA cAMPi + ACA 


kgXi 


1^5 9 = +1 


kg = 0.306min"^ 


RegA + cAMPi ^ RegA 


kloX4X5 


fS 10 = -1 


fcio = 0.816^iM" Vin"^ 


ACA cAMPe + ACA 


kiixi 


1^6 11 = -1-1 


fell = 0.686min"^ 


cAMPe 


kl2X6 


1^6 12 = —1 


ki2 = 4.998min"^ 


cAMPe ^ CARl + cAMPe 


kvjxe 


1^7 13 = +1 


ki3 = 22.54min"^ 


CARl ^ 


k\4,X7 


U7 14 = —1 


ki4 = 4.59min"^ 



TABLE I. The reactions and their associated rate constants, given in micromoles and minutes. All the reactions have mass- 
action kinetics. Only non-zero entries of the stoichiometry matrix, Ui^, are shown. 
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Some results from the one-cell model will be shown, before looking at the two cell model. The reaction pa- 
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P{cj) Im(C56(w)) 




Re(C56(a;)) 



FIG. 4. Theoretical power spectra for the internal (larger 
peak) and external cAMP in the one-cell model. The peak is 
around u} = 0.83. The reaction rates are given in Table 1. In 
[18| the cell volume was chosen to be 3.672 x 10""l. 



FIG. 6. Parametric plot of the (theoretical) CCF for the in- 
ternal and external cAMP fluctuations in the one-cell model. 
Only results for 0.5 < to < 1.2 are shown. The reaction rates 
are given in Table 1. 
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FIG. 5. Magnitude of the (theoretical) CCF for the internal 
and external cAMP fluctuations in the one-cell model. The 
reaction rates are given in Table 1. 



cAMPi 
20 000 h 




FIG. 7. Time series for cAMPi (in particle numbers) in differ- 
ent cells for a three-cell model, generated using the Gillespie 
algorithm. The three species are given different initial con- 
ditions, but within a short amount of time they oscillate in 
phase with each other around their common fixed point. 



rameters are chosen to be those which yielded the steady 
state in Figure IB in 'TS]. As in Section 2, we use the 
linear noise approximation to quantify the fluctuations 
around the steady state. We can then use Eq. ([t]) to cal- 
culate the PSDM and Eq. ^ to calculate the CCF. The 
specific form of the analytic results are not very illumi- 
nating and we content ourselves with showing the results 
graphically. Figure |4] shows the power spectra for the 
internal and external cAMP. The values of the reaction 
parameters are given in the caption. Figure [5] shows the 
magnitude of the CCF for these two species, showing a 
very strong correlation, especially in the frequency range 
within which the oscillations are significant. The CCF 
is displayed parametrically in Figure |6] It has both real 
and imaginary components, which indicates a phase lag. 

In [18], the authors construct models of many cells, 
which are all coupled via the external cAMP. A diagram 
of such a model, with three cells, is shown in Figure 5A 



of 18 . If the reaction parameters within each cell are 
identical, the oscillations in the three cells become syn- 
chronised with zero-phase lag. This effect can be easily 
seen by looking at the time series. This is shown for a 
three-cell model in Figure [Tj To apply our analysis to 
a multi-cell model we proceed as we did for the one cell 
model: for an n cell model there will be 6n -f 1 variables. 

Synchronisation with a non-zero phase lag can be 
found when the reaction parameters differ from cell to 
cell. We illustrate this with a two cell model, where the 
reaction parameters in one cell are those used for the 
single cell model, and the reaction parameters in the sec- 
ond cell are obtained by making random perturbations 
to the original parameters. Below the relation between 
the internal cAMP oscillations in each cell is examined, 
for one such parameter choice. Figure plshows the power 



spectra for these two species. Figures[9^ & 10 show the 
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P(co) 




0.5 1.0 1.5 2.0 

FIG. 8. Theoretical power spectra for cAMPi in each cell of 
the two-cell model. The dashed line corresponds to the oscil- 
lations in the second cell, for which the reaction parameters 
are: fci = 2.3min~^, k2 = 0.96/iM~'^min~^, kz = 2.1min~^, 
ki = 1.9min"\ fcs = 0.4min"\ ke = 0.72^iM" Vin"\ 
k-j = 0.7AiMmin-\ fcg = LOS/iM-^min-^ fcg = 0.26min-\ 
fcio = 0.89^iM" Vin"\ fcn = 0.46min~\ fcia = 15min"\ 
fci4 = 5.8min^^. The parameter fci2 was set to 9min^^. 
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FIG. 10. Parametric plot of the CCF for the cAMPi in each 
cell of the two-cell model. The solid line is the theoretical 
result, the dots are from 6000 numerical simulations. Results 
are shown for < a; < 3. 
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FIG. 11. Phase lag for oscillations of cAMPi in each cell for 
the two-cell model. The solid line is the theoretical result, the 
dots are from 6000 numerical simulations. Over the frequency 
range where oscillations are significant, the phase lag remains 
fairly constant. 



FIG. 9. The magnitude of the CCF for the cAMPi in each 
cell of the two-cell model. The solid line is the theoretical 
result, the dots are from 6000 numerical simulations 



relevant CCF for these two species, whilst Figure [TT dis- 
plays the phase lag present, calculated using Eq. (10 1. 
The theoretical results show good agreement with those 



from numerical simulation. Although the reaction pa- 
rameters are different in each cell, there remains a strong 
shared signal in the oscillations. 

Throughout this section so far, we have assumed that 
all of the cells have the same volume. We will now look at 
the consequences of relaxing this restriction, which could 
be due to natural variation in the cell size. We find that 
a phase lag is introduced when the cells are of different 
sizes. Given that, for this particular system, it is de- 
sirable to have in phase synchronisation, it is important 
to see how large these phase lags are for various differ- 



ences in cell size. To do this, we again studied a two-cell 
model, in which we fixed the volume of one cell, at the 
volume used previously which we will call Vi , and the ex- 
tracellular region (fixed at 5V7), and varied the volume of 
the other cell. The reaction parameters in each cell were 
identical and were chosen to be those used for the one 
cell model. The parameter governing the rate of degra- 



dation of cAMPe, ki2, was set to 12min~ . Figure [12 
shows details of the CCF for the case where the second 
cell has volume lAVj. The parametric plot shows that 
the imaginary component of the CCF is much smaller 
than the real component: this means that the phase lag 
is very small. We then increased the volume of the sec- 
ond cell to 2Vi. Figure [l3| shows that the imaginary part 
of the CCF is larger than before, indicating a more sig- 
nificant phase lag. The phase spectrum for this system 
is shown in Figure |14[ At the frequency for which the 
oscillations are most significant (the frequency for which 
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FIG. 12. Results for a two-cell model, where the cell volumes 
are not identical. Here, the cell volumes were chosen to be 
Vi and lAVi, where Vi is the cell volume chosen by Kim. et. 
ai. The reaction parameters were the same in each cell, and 
were those used for the one-cell model. Top: The magnitude 
of the CCF for the cAMP fluctuations in each cell. Bottom: 
A parametric plot, showing the real and complex parts of the 
same CCF for the frequency range < oj < 2. A small phase 
lag has appeared, due to the different cell volumes. 
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FIG. 13. Results for a two-cell model, where the cell volumes 
are not identical. Here, the cell volumes were chosen to be Vi 
and 2Vi, where Vi is the cell volume chosen by Kim. et. al. 
The reaction parameters were the same in each cell, and were 
those used for the one-cell model. Top: The magnitude of 
the CCF for the cAMP fluctuations in each cell. Bottom: A 
parametric plot, showing the real and complex parts of the 
same CCF for the frequency range < w < 2. A phase lag 
has appeared, due to the different cell volumes. 



the povifer is greatest), the phase lag is about 0.2 radians, 
roughly 11°. This shows that quite large differences in 
volume are required to introduce measurable phase lags: 
varying the cell size by a few percent does not have a sig- 
nificant effect. This is positive from the point of view of 
in-phase synchronisation being robust in this system. We 
also repeated this analysis whilst varying the volume of 
the extracellular region. However, this did not produce 
significantly different results from those presented here. 

In this section we have studied a model of cAMP os- 
cillations in Dictyostelium, due to Kim et. al. Using the 
theoretical framework of the LNA, we have shown how 
analytical results can complement the numerical results 
which have already been obtained from this model TJJ. 
The analysis performed in this section extends naturally 
to models containing many cells. However, in reality 
the Dictyostelium cells can form aggregates containing 
thousands of cells. In this case, the distances between 
cells will be much larger, and the diffusion speed of the 
molecules would need to be taken into consideration. A 
spatial model would be required to capture this effect, 
although the mechanism by which the cells synchronise 
would be the same as the one outlined here. 




FIG. 14. The phase spectrum for the cAMP fluctuations 
in each cell in a two-cell model for the case when the cell 
volumes are not identical. The cell volumes were Vi and 2Vi. 
The power spectra for these oscillations peak at c^; = 0.9. The 
phase lag for oscillations of this frequency is 0.2 radians. 



IV. DISCUSSION 

In this work we have shown how stochastic oscillations 
in biochemical models can synchronise, and how analyt- 
ical expressions for this effect may be found. Compar- 
isons between numerical results and theoretical predic- 
tions for the CCF (and related quantities) were found to 
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be good, especially in frequency ranges where the oscil- 
lations were significant. There have been very few mod- 
els devised which study the synchronisation of stochastic 
oscillations across multiple cells, and those which have 
been constructed, such as [18j , are examined using nu- 
merical simulation. We hope that the analytical work 
described here can complement the numerical approach. 
In our study of the model of cAMP signalling, we found 
that the stochastic oscillations across cells strongly influ- 
enced each other, and oscillations of species in different 
cells synchronised rapidly, as illustrated in Figure [7] In 
the case where the cells were of equal volume, and the 
reaction parameters across cells were the same, the oscil- 
lations synchronised in phase. Introducing changes to the 
reaction parameters in different cells, or making the cell 
volumes different to each other introduced a phase lag. 
These findings are similar to those reported in [22], in 
the context of epidemiological modelling. However, quite 
large changes had to be made before the lag was found 
to be significant. These results are positive for in-phase 
synchronisation being robust in this system. One slightly 
surprising result we found was that the size of the phase- 
lag was unaffected by varying the volume of the extra- 
cellular compartment. This is probably not realistic, and 
could be due to the assumption that all compartments 
are well mixed. For a larger extracellular compartment, 
communication between cells would be expected to take 
longer, which could affect the phase lag. 

The analysis used to study the model of cAMP os- 
cillations was also applied to a model of oscillations in 



yeast glycolysis, due to Wolf & Heinrich pTj. Although 
not presented here, very different results were obtained 
for this model. The magnitude of the CCF for oscillators 
in different cells was extremely small. Using the reaction 
parameters chosen in the original article, we obtained 
CCFs with a magnitude of < 0.01 over the relevant fre- 
quency range. This indicates that the shared signal in 
the oscillators is extremely weak. The magnitude of the 
CCF remained small throughout an extensive parameter 
search of the model. This weak coupling between the 
cells could explain why, in the deterministic framework, 
the oscillators take a longer-than-expected time to syn- 
chronise. The relation between the time to synchronise 
and the magnitude of the CCF could be an interesting 
avenue for further work. The magnitude of the fluctu- 
ations, here governed by the system-size, could also in- 
fluence whether synchronisation is observed i.e. if the 
system-size is small (and therefore the fluctuations are 
relatively large), is the synchronisation of two oscillators 
observed if the magnitude of the corresponding CCF has 
low or intermediate value? Such questions we leave for 
the future. Our aim here has been to formulate the ana- 
lytic study of synchronisation of stochastic oscillators in a 
biochemical context, and to illustrate its use in a specific 
model. We hope that this will encourage others use these 
techniques to investigate other biochemical systems. 
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